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Abstract — A  recent  study  [1]  suggests  that,  under  certain 
circumstances,  concordant  changes  in  cerebral  intravascu¬ 
lar  oxygenation  and  mean  arterial  blood  pressure  reflect 
impaired  cerebrovascular  autoregulation.  In  this  paper, 
we  propose  a  new  measure  to  quantitate  this  concordance, 
derived  from  the  common  subspace  of  these  signals.  The 
method  is  compared  to  correlation  and  coherence  analysis 
with  respect  to  our  application,  but  it  is  also  suited  for 
other  biomedical  signals.  Furthermore,  this  model-based 
approach  is  not  restricted  to  applications  involving  only  2 
signals. 

Keywords — subspace-based  modeling,  correlation,  coher¬ 
ence,  cerebrovascular  autoregulation,  near-infrared  spec¬ 
troscopy 


NIRS,  and  MAP,  assessed  by  intravascular  catheterisa¬ 
tion.  If,  in  this  way,  infants  at  high  risk  for  GMH- 
IVH/PVL  can  be  identified  before  the  occurence  of  the 
lesions,  it  might  be  possible  to  correct  the  cerebral  circu¬ 
latory  disturbance  and  prevent  the  lesions. 

The  mathematical  problem  is  thus  the  quantitation  of 
the  concordance  between  HbD  and  MAP.  For  this  pur¬ 
pose,  three  different  measures  are  described  and  com¬ 
pared.  Furthermore,  some  preliminary  results  obtained 
from  measured  data  are  given. 

II.  Methods 


I.  Introduction 

The  most  important  forms  of  brain  injury  in  pre¬ 
mature  infants  (germinal  matrix-intraventricular  hem¬ 
orrhage  (GMH-IVH)  and  periventricular  leukomalacia 
(PVL))  are  caused  in  considerable  part  by  disturbances  in 
cerebral  blood  flow  (CBF).  Premature  infants  have  a  par¬ 
ticular  propensity  for  development  of  disturbances  in  CBF 
for  2  major  reasons.  First,  because  alterations  in  mean 
arterial  blood  pressure  (MAP)  are  very  common  in  such 
infants.  And  second,  because  cerebral  autoregulation,  the 
mechanism  by  which  CBF  is  maintained  constant  despite 
alterations  in  MAP,  is  either  defective  or  absent  in  at  least 
some  infants.  So,  there  is  a  relationship  between  impaired 
autoregulation  and  the  occurence  of  brain  injury. 

Two  studies  [2]  [3]  showed  a  relationship  between 
CBF  and  MAP  and  concluded  that  autoregulation  may 
be  defective  in  some  premature  infants.  Different  au¬ 
thors  did  not  find  any  correlation  between  CBF  and 
MAP,  however  they  used  singular  measurements  and  not 
continuous  measurements.  A  recent  study  suggests  [1] 
that  concordant  changes  in  cerebral  intravascular  oxy¬ 
genation  (HbD)  and  MAP  reflect  impaired  cerebrovas¬ 
cular  autoregulation  in  continuous  measurements.  A 
previous  study  [4]  showed  a  strong  correlation  between 
HbD,  measured  non-invasively  by  near-infrared  spec¬ 
troscopy  (NIRS)  as  the  difference  between  the  concen¬ 
tration  changes  of  oxygenated  hemoglobin  (L/hC^)  and 
deoxygenated  hemoglobin  ( Hb ),  and  volemic  CBF,  deter¬ 
mined  by  radioactive  microspheres  if  the  arterial  oxygen 
saturation  (SaO 2)  does  not  change  appreciably  during  the 
measurement.  Thus,  under  such  conditions,  HbD  is  a 
measure  of  CBF. 

Consequently,  premature  infants  with  impaired  cere¬ 
brovascular  autoregulation  could  be  identified  by  simul¬ 
taneous,  continuous  measurements  of  HbD,  assessed  by 


A.  Correlation 


The  most  straightforward  way  to  quantitate  the  concor¬ 
dance  between  two  signals  x(t)  and  y(t)  is  the  correlation 
coefficient  COR ,  estimated  as: 


COR  = 


Z?=Mt)-x){y(t)-y) 

y/Et=Mt)-^)2Et=Mt)-y)2 


(i) 


where  N  is  the  number  of  samples  considered  and  x ,  y 
stands  for  the  mean  of  x(t),  respectively  y(t).  If  one  is 
only  interested  in  the  correlation  in  a  specific  frequency 
band,  the  signals  should  be  filtered  with  a  band-pass  filter 
containing  the  frequency  band  of  interest  before  calculat¬ 
ing  the  correlation  coefficient. 


B.  Coherence 


Another  way  to  quantitate  the  correlation  in  a  fre¬ 
quency  specific  manner  is  based  on  the  coherence  func¬ 
tion,  defined  as: 


cxy{f) 


\P*y(f)\2 

Pxx{f)Pyy{f ) 


(2) 


where  Pxy{f )  is  the  cross  power  spectral  density  of  x(t) 
and  y(t),  and  Pxx{f),  Pyy{f )  the  power  spectral  densities 
of  x(t),  respectively  y(t).  The  spectral  density  functions 
are  estimated  using  Welch’s  method  [5]  as  follows: 

•  x(t),  y(t)  are  divided  in  overlapping  sections  of  n  points; 

•  each  section  is  detrended  and  windowed; 

•  the  length  n  fast  Fourier  transform  (FFT)  is  calculated 
of  each  section  (these  FFTs  are  called  periodograms); 

•  the  pointwise  products  of  the  spectra  of  x(t)  and  y(t) 
are  averaged  over  the  overlapping  sections  to  form  Pxy(f); 

•  the  pointwise  squared  FFTs  of  x(t)  and  y(t)  are  aver¬ 
aged  over  the  sections  to  form  Pxx{f ),  respectively  Pyy{f)\ 

•  Cxy  is  calculated  by  means  of  (2). 
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The  coherence  function  is  a  measure  of  the  linear  depen- 
dance  between  the  signals  x{t)  and  y(t)  at  a  given  fre¬ 
quency.  A  coherence  of  1  indicates  perfect  frequency  spe¬ 
cific  correlation,  a  coherence  0  indicates  a  complete  lack  of 
frequency  specific  correlation.  To  obtain  a  measure  of  the 
coherence  over  a  specific  frequency  band,  the  coherence 
function  is  averaged  over  that  frequency  band. 

C.  HTLS-SEP 


Then,  to  determine  the  common  subspace  C  of  Hx  and 
Hy ,  matrices  X  and  Y  must  be  calculated  such  that: 


HxX  =  HyY  =  C 


Since  this  condition  can  also  be  written  as: 


[Hx  -  Hy] 


X 

Y 


=  0 


(8) 

(9) 


This  method  was  first  described  in  [6]  and  is  briefly 
described  in  this  section.  The  goal  of  the  HTLS-SEP 
method  is  to  fit  one  linear  model  to  two  signals.  Only  if 
the  two  signals  have  the  same  dynamical  properties  (or,  in 
other  words,  have  common  signal  poles),  the  determined 
system  will  be  able  to  explain  both  signals  sufficiently 
well.  Let  us  start  from  the  hypothesis  that  the  two  sig¬ 
nals  have  the  same  dynamical  properties.  Mathematically 
this  is  reflected  in  the  proposed  model  of  x(t)  and  y(t)  as 
follows: 

K 

x(t)  =  ^2  clz{  +  m(i),  t  =  0, 1, . . . ,  N  -  1  (3) 

k= 1 

K 

y(t )  =  clzl  +n2(t),  t  =  0,l,...,N  -1  (4) 

fc= i 

where  ni(t)  and  n2(t)  represent  the  noise,  and 
Zk,  k  =  1, . . . ,  K  the  so-called  signal  poles: 

Zk  =  eU2nfk-dk)At  (5) 

where  fk,  k  =  1  represent  the  frequencies, 

dk,  k=l, . . . , K  the  dampings  and  At  the  constant  sam¬ 
pling  interval;  furthermore,  ck,  k  =  1, K  represent 
the  so-called  complex  amplitudes  associated  with  x(t): 


it  is  clear  that  X  and  Y  can  easily  be  derived  from  the 
null  space  of  the  block  Hankel  matrix  [Hx  Hy\ . 

The  common  poles  are  thus  estimated  from  the  common 
subspace  C  by  means  of  the  HTLS  algorithm1  [8],  giving 
the  pole  estimates  zk,  k  =  1  This  is  possible 

because  the  column  spaces  of  Hx  and  Hy  possess  the  shift 
invariance  property.  Multiplying  them  with  a  matrix  to 
the  right  does  not  change  this  property. 

After  estimation  of  the  common  poles,  the  phases  and 
the  amplitudes  are  estimated  for  x(t)  and  y(t)  separately 
as  the  least  squares  solutions  to  (3),  respectively  (4),  with 
zk  replaced  by  the  estimates  zk. 

Using  (3)  and  (4)  with  all  parameters  replaced  by  their 
estimates,  we  reconstruct  the  part  of  x(t)  common  to  x(t) 
and  y(t),  giving  xc{t).  The  residual  xr(t)  can  then  be 
defined  as: 

xr(t)  =  x(t)  —  xc(t)  (10) 

Note  that  if  x(t )  and  y(t)  do  not  have  many  components  in 
common,  xc(t)  will  be  a  very  bad  fit  to  x(t).  Therefore, 
the  ratio  of  the  energy  of  the  common  part  to  the  sum 
of  the  energy  of  the  common  part  and  the  energy  of  the 
residual,  is  used  as  a  measure  of  the  importance  of  the 
common  part  in  the  original  signal  x(t): 

xr(t)  =  x(t)  —  xc(t)  (11) 


c%  =  a%e (6) 

where  ak,  k  =  1  stands  for  the  amplitude  and 

< f>k ,  k  =  1, . . . ,  K  for  the  phase.  In  a  similar  way  the  com¬ 
plex  amplitudes  ck,  k  =  of  y(t)  are  defined  as 

a  function  of  avk,  k  =  1, . . . ,  K  and  c t>vk ,  k  =  1, . . . ,  K. 

From  (3)  and  (4),  it  is  clear  that  the  two  signals 
x(t)  and  y(t)  are  modeled  with  the  same  signal  poles 
zk,  k  =  1,  2, . . . ,  K  (meaning  that  they  share  the  dynam¬ 
ical  properties  of  one  linear  system)  but  with  different 
complex  amplitudes.  To  determine  the  common  poles, 
a  subspace-based  approach  is  used.  First,  Hankel  matri¬ 
ces  Hx  and  Hy  are  constructed  from  x(t),  respectively 
y(t),  as  follows: 


CPCX 


— 1 
2^t= o 


(®c(*))2  +  Et=o1 


(xr(tW 


(12) 


If  the  signals  consist  only  of  common  poles  than  CPCX 
approaches  one,  if  there  are  no  common  poles  CPCX  will 
be  close  to  zero.  If  the  signals  are  normalized  (standard 
deviation=l),  CPCX  can  also  be  regarded  as  a  measure  of 
the  mean  square  error  of  xc(t)  with  respect  to  x(t),  since: 


CPCX  «  1  (13) 

t= 0 

In  a  similar  way  yc(t ),  yr(t)  and  CPCy  can  be  defined. 


x(o) 

x{l)  . 

..  x(N—m+ 1) 

x(l) 

x(2)  . 

. .  x(N—m+ 2) 

x(m—  1)  x(m) 

x(N-l) 

^^HTLS  is  a  subspace-based  pole  estimation  method  that  exploits 
the  shift  invariance  property.  It  is  an  alternative  to  the  TLS- 
ESPRIT  method  [7]:  HTLS  uses  the  singular  value  decomposition 
of  the  data  matrix,  whereas  TLS-ESPRIT  uses  the  eigenvalue  de¬ 
composition  of  the  sample  covariance  matrix. 
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III.  Experimental  results 

A.  Measurements 

From  several  premature  infants,  MAP,  SaO 2  and  HbD 
were  measured  simultaneously  at  the  University  Hospital 
Leuven.  HbD  was  calculated  as  the  difference  between  the 
concentration  changes  of  U6O2  and  Hb.  These  changes 
were  acquired  using  near-infrared  spectroscopy  (NIRO- 
300®,  Hamamatsu).  The  data  are  recorded  at  a  sampling 
frequency  of  100  Hz  by  a  data  acquisition  system  Codas 
(CODAS®,  Dataq  Instruments,  USA)  and  stored  on  a 
PC.  The  MAP-  and  SaC>2-  signals  are  analog  and  digitized 
afterwards  by  the  CODAS-system.  The  NIRO-300  signals 
are  digital  with  a  sampling  rate  of  6  Hz.  They  are  con¬ 
verted  to  analog  signals  with  a  sample-and-hold  function 
before  their  introduction  in  the  CODAS-system.  From 
this  100  Hz  data,  the  average  values  for  non-overlapping 
5-second  intervals  were  calculated  (0.2  Hz).  This  sampling 
frequency  is  still  high  enough,  since  major  physiological 
importance  can  be  attributed  to  the  frequency  band  of 
0-0.01Hz  (i.e.  changes  occuring  over  several  minutes)  [1]. 

B.  Results 

Since  the  concordance  between  the  signals  might  vary 
as  a  function  of  time,  a  sliding  window  approach  was 
used:  the  presented  measures  were  calculated  on  30- 
minute  recordings  (i.e.  360  samples),  and  plotted  against 
the  starting  time  of  the  30-minute  window.  In  this  way  it 
is  possible  to  track  the  changes  of  the  measures  over  time. 

For  the  correlation  coefficient,  we  used  the  5-second  av¬ 
eraged  values  without  filtering.  When  the  signals  are  low- 
pass  filtered  to  the  frequency  range  of  interest  (0-0.01Hz), 
the  correlation  coefficient  is  only  slightly  higher  in  highly 
correlated  windows,  but  also  higher  in  uncorrelatecl  win¬ 
dows.  Therefore,  the  correlation  coefficient  of  the  un¬ 
filtered  signals  is  more  sensitive  to  discriminate  between 
correlated  and  uncorrelated  windows. 

Several  parameters  have  an  influence  on  the  estimation 
of  the  spectra  for  the  calculation  of  the  coherence  func¬ 
tion.  In  the  first  place,  the  length  n  of  the  sections,  used 
for  the  calculation  of  the  FFTs,  and  the  number  of  over¬ 
lapping  points  between  consecutive  sections.  The  smaller 
the  value  of  n,  the  more  periodograms  can  be  averaged 
and  hence  the  lower  the  variance  of  the  estimates  of  the 
spectra.  On  the  other  hand,  the  smaller  n,  the  lower  the 
frequency  resolution.  Experiments  with  varying  n  showed 
that  n  =  144  gives  the  best  results  (the  frequency  resolu¬ 
tion  is  then  0. 0014Hz).  In  order  to  average  as  much  as 
possible,  we  used  maximum  overlap  between  the  sections 
(i.e.  n—  1  =  143  data  points).  The  average  of  the  coherence 
function  over  the  frequency  band  of  interest  (0-0.01Hz) 
was  used  as  measure,  and  is  called  COH  in  the  remainder 
of  this  paper.  Before  calculating  the  FFT  of  each  section, 
the  mean  was  removed  and  a  Hanning  window  applied. 

The  main  problem  for  the  HTLS-SEP  method  is  the 
determination  of  the  model  order  of  the  common  sub¬ 
space  K,  which  is  twice  the  number  of  frequency  com¬ 


ponents.  Since  K  is  not  known  a  priori,  the  following 
approach  is  used.  First,  x(t )  and  y(t)  are  modeled  with 
I\=2.  Then,  as  long  as  ex,  being  the  mean  of  |a;(f)  —  a:(t)| 
over  the  30-minute  window,  and  ey ,  being  the  mean  of 
|y(f)  —  y(t)  |  over  the  same  period,  decreases,  K  is  in¬ 
creased  with  2.  If  either  ex  or  ey  increases,  the  algorithm 
is  stopped.  The  model  of  HbD  with  the  previous  value  of 
K  is  then  used  to  calculate  CPC  by  means  of  (12).  Be¬ 
fore  modeling,  the  signals  of  the  30-minute  windows  were 
normalized  (such  that  mean=0,  standard  deviation=l). 

An  example  of  the  measurements  on  a  premature  in¬ 
fant  are  shown  in  Fig.  1,  together  with  the  corresponding 
measures.  The  total  length  of  the  measurement  was  four 
hours.  The  MAP  signal  shows  only  modest  fluctuations 
during  the  first  two  hours.  There  is  also  modest  variabil¬ 
ity  in  HbD  during  the  first  two  hours,  except  between 
the  40th  and  60th  minute,  where  HbD  clearly  decreases. 
However,  this  decrease  is  accompanied  by  a  decrease  in 
SaO 2,  as  can  be  seen  in  Fig.  2.  Fig.  2  also  shows  MAP 
and  HbD  during  the  30-minute  window  containing  this  de¬ 
crease,  together  with  the  models  of  MAP  and  HbD.  Since 
the  condition  of  constant  Sa02  does  not  hold,  HbD  can 
not  be  assumed  to  be  proportional  to  CBF  during  this  pe¬ 
riod,  and  therefore  the  results  should  be  interpreted  care¬ 
fully.  The  analysis  of  this  30-minute  window  results  in  low 
COR  and  COH  values  (0.26  and  0.39  respectively);  on  the 
other  hand,  CPC  is  high  (0.86).  The  latter  indicates  that 
the  dynamics  of  the  Sa02  is  not  only  reflected  in  HbD, 
but  also  in  MAP  although  to  a  lesser  degree.  A  possi¬ 
ble  explanation  of  this  finding  is  that  the  patient  might 
have  had  an  open  ductus  Botalli.  This  causes  decreases 
in  oxygenation  which  are  often  accompanied  by  decreases 
in  blood  pressure.  However,  further  research  is  necessary 
to  clarify  this  hypothesis.  During  the  last  two  hours  of 
the  recording,  large  (spontaneous)  fluctuations  in  MAP 
are  associated  with  parallel  changes  in  HbD,  while  Sa02 
remains  constant.  This  is  a  case  of  impaired  cerebrovas¬ 
cular  autoregulation.  During  this  period,  CPC  (close  to 
1)  and  COR  (about  0.8)  are  higher  than  COH,  suggesting 
that  CPC  and  COR  are  more  sensitive  in  order  to  detect 
impaired  autoregulation.  To  evaluate  the  predictive  value 
of  the  measures  with  respect  to  subsequent  brain  injuries 
statistically,  a  larger  follow-up  study  is  necessary. 

IV.  Conclusions 

A  new  measure  to  quantitate  how  similar  the  dynamics 
in  two  (or  more)  signals  are,  was  presented.  It  is  derived 
from  the  common  subspace  of  the  signals.  The  measure 
was  used  to  quantitate  the  concordance  between  cerebral 
intravascular  oxygenation  and  mean  arterial  blood  pres¬ 
sure,  which  reflects  impaired  cerebrovascular  autoregula¬ 
tion.  The  recordings  analyzed  so  far  indicate  that  the  new 
method  and  correlation  are  better  measures  to  detect  im¬ 
paired  autoregulation  than  coherence  analysis.  A  larger 
follow-up  study  is  needed  to  confirm  the  performance  of 
the  new  measure,  and  to  evaluate  the  predictive  value 
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Fig.  1.  Simultaneous  changes  in  MAP  (a)  and  HbD  (b)  in  a  premature  infant,  with  corresponding  measures  (c):  COR  (dashed  line),  COH 
(dash-dotted  line)  and  CPC  (solid  line). 


of  the  different  measures  with  respect  to  possible  brain 
damage  statistically.  Furthermore,  the  new  measure  truly 
looks  for  similarity  of  the  dynamics  between  signals  and 
should  therefore  be  better  suited  to  see  whether  the  sig¬ 
nals  are  generated  by  one  and  the  same  biomedical  pro¬ 
cess.  Another  advantage  of  the  new  method  is  that  it 
allows  to  analyze  the  dynamical  properties  of  more  than 
two  signals  simultaneously.  This  will  allow  us  in  the  future 
to  include  other  relevant  physiological  parameters  like  the 
arterial  carbon  dioxide  tension  ( PCO2 )  in  the  analysis. 
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